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ABSTRACT 

In this report, we have designed an essentially non-oscillatory reconstruction for functions 
defined on finite-element type meshes. Two related problems are studied : the interpolation 
of possibly unsmooth multivariate functions on arbitrary meshes and the reconstruction 
of a function from its average in the control volumes surrounding the nodes of the mesh. 
Concerning the first problem, we have studied the behaviour of the highest coefficients of the 
Lagrange interpolation function which may admit discontinuities of locally regular curves. 
This enables us to choose the best stencil for the interpolation. The choice of the smallest 
possible number of stencils is addressed. Concerning the reconstruction problem, because of 
the very nature of the mesh, the only method that may work is the so called reconstruction 
via deconvolution method. Unfortunately, it is well suited only for regular meshes as we 
show, but we also show how to overcome this difficulty. The global method has the expected 
order of accuracy but is conservative up to a high order quadrature formula only. 

Some numerical examples are given which demonstrate the efficiency of the method. 


*This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-18605 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1 Introduction 


During the past few years, a growing interest has emerged for building high order accurate 
schemes (i.e of order greater that 2) for compressible flows simulations. It is well known that 
even for smooth initial conditions, these flows may develop discontinuities that make linear 
schemes useless. 

At the beginning of the 80’s, the class of Total Variation Diminishing schemes appeared 
and they have been successfully and widely used with many types of meshes (see for example, 
[1] for a review and, among many others, [2] for simulations on finite element type meshes). 
Nevertheless, one of their main weaknesses is that the order of accuracy falls to first order 
in regions of discontinuity and at extrema, leading to excessive numerical dissipation. 

Various methods have been proposed to overcome this difficulty (adaptation of the mesh 
for example) but one promising way may also be the class of the Essentially Non-Oscillatory 
schemes (E.N.O. for short) introduced by Harten, Osher and others [3, 4, 5, 6, 7]. 

The basic idea of E.N.O schemes is to use a Lagrange type interpolation with an adapted 
stencil : when a discontinuity is detected, the procedure looks for the region around this 
discontinuity where the function is the smoothest. Then a reconstruction technique may be 
applied which enables approximation of the function to any desired order of accuracy from 
its average in control volumes surrounding the mesh points. The approximation is done so 
that it is conservative. 

Some attempts have been made to extend these ideas to multidimensional flows (see for 
example [8]), but only for structured meshes. 

In this report, we intend to study the problem of the reconstruction, up to any order 
of accuracy, of a given function given either by its value at the nodes of a triangulation or 
by its averages on control volumes defined around these nodes so that, in the second case, 
the reconstruction is conservative. This latter problem has already been studied, for smooth 
functions only, by Barth et al. [9] but their method does not appear to generalize easily to 
unsmooth functions. 

The outline of this report is as follows. In the first part, we give some basic facts 
about Lagrange interpolation in several dimensions. In particular, we study the problem of 
the localization of regions of smoothness from the Lagrange interpolation coefficients that 
generalize those known in one dimension. These results seem to be new. Then, we propose 
an algorithm for E.N.O. interpolation and we try to give some indications for selection of the 
smallest family of the possible stencils for second and third order approximation, the only 
ones considered in this report. We also propose an adaptation of the so-called reconstruction 
via deconvolution procedure that was originally built for regular meshes, and indicate why, 
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in general, the conservation property must be lost to ensure a high order of approximation. 
Some numerical tests indicate the performance of this method. 

1.1 Notations 

• IRnpi, Y] : finite dimensional vector space of two variable polynomials over IR, 

• N(n) = ( n — ^ : dimension of IR n [A', Y], 

• admissible stencil for solving the Lagrange problem in IR„[X, Y], see section 2.1.2, 

• | It/’ll is the euclidian norm of U, 

• Diju = @ - where i + j- l, 

3 a xo 3 y 

• D n u : n-th derivative of u. 

2 Lagrange interpolation on arbitrary sets of IR 2 

2.1 Sets of admissible points 

2.1.1 The polynomials in IR 2 

We will denote by IR[A, Y] the vector space of the polynomials of two variables ( X and Y) 
with coefficients in IR. An element of IR[X, Y] may be described by its (finite) expansion in 
terms of powers of X and Y : 

P(X,Y) = J2 £ a, jJTYi (1) 

Z=1 i+j=l i,j> 0 

The highest integer such that at least one of the coefficients of the monomials X'Y 3 is non 
zero is called the total degree of P. 

If (xo,yo) is a point of IR 2 , another expansion of P may be written in terms of the 
monomials ( X — x o ) J ( Y — Vo) 3 with the help of the Taylor formula. The total degree of P 
does not depend on the point (a; 0 ,yo)- 

In the sequel, we will denote by JRn[X,y] the (finite) vector space of the polynomials 
of IR[X, Y] with total degree less or equal to n. This vector space has dimension N(n) = 
( w — a basis of which is the set of monomials (X - x 0 )'(Y — y 0 ) 3 of total degree 

i + j less or equal to n. 

Let us now describe another interesting basis of 1R„[.A, Y]. Consider (A, B , C ) a triangle of 
IR 2 and let us denote by A^, Ag, Ac the barycentric coordinates of the three points (A, B, C ) 
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defined, for any point M, by : 


M — A + Ag B 4 * A c C 
A a + A# + Ac = 1 


( 2 ) 


It is easy to se that, for any pair of points, say A and B, the set |A^ A^| 
basis of IRnpf, Y]. 


*+j<n 


is also a 


2.1.2 The Lagrange interpolation problem in IR 2 

The Lagrange interpolation problem may be formulated as follows : 


Given N and n, two integers, a family of N points in IR 2 , .S^") = (Ai) 1<i<N and 
N real values find an element P o/ILJX, Y] such that for any point 

Ai ofSM, one has P(Ai ) = tq. 


In the sequel, we will often make no distinction between an element of S^ n \ Ai, an d its 
coordinates in a suitable frame, ( x\,yi ). 

For this problem to have a solution, two conditions must be fulfilled : 

1. one must have N = ^ 

2. the following generalized Van der Monde determinant must be non zero : 


A S (n) 


det x) yj | 


i + j < n 


1 Xi t /1 x" x" 1 yi •••x? 

1 x N y N ••• xft x^ -1 j/7v ”’Xn 


( 3 ) 


We will say that the set is admissible if A 5 (n) ^ 0. If a set is admissible, there 
exist ^ ~ ^ coefficients , (a, j), such that the solution of the Lagrange problem is : 


p = E E ». (4) 

1=1 i+j=l i,j>0 

The problem of characterizing the admissible sets has been widely studied, see [10] for 
example and the references therein. 


Remarks : 

1. The condition (3) has been given for the basis X'Y J of IR^fX, Y]. A similar and 
equivalent condition could have been given for the two other bases we have mentioned. 
The formula (4), provided that the monomials X'Y* are replaced by the elements of 
the new basis, is also true. 
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2. If card <SH) = 3, this condition is nothing more than the one which says that the three 
points must not be aligned. 

3. If we were in IR, this determinant would be the classical van der Monde determinant. 

4. The set of ~ -uplets where the condition (3) is not fulfilled is an algebraic 

curve of 1R/ V ^ and consequently a closed subset of measure zero in IR^'”'. 

In the next section, we address the question of the practical calculation of the coefficients 


2.1.3 Determination of the Lagrange expansion 

In this section, we use the monomials X'Y j for expanding polynomials, but any other basis 
would be suitable and the results are immediately transferable. 

The coefficients of the Lagrange expansion are the solution of the linear N(n ) x N(n ) 
system : 

w/ = 53 13 a * yi for a11 ( x i>vO G £ (n) - ( 5 ) 

1=1 i+j=l i t j > 0 

Since condition (3) is true, the Cramer formula applied to (5) gives the answer. 

Several authors have tried to generalize the Newton formula that make the Lagrange 
interpolation efficient from a numerical point of view, and a very general answer has been 
given by Muhlbach [11, 12]. 

In these papers, he addresses the problem of the Lagrange interpolation by a set of 
functions (fi) ieJ on a set of points (A,-) igJ . He calls the set (/,) a Cebyscv-system if given 
any function /, for any pair of subsets of 7, L and M, having the same (finite) number of 
elements, there exist real numbers a/ such that : 


«.■ = 53 OTj/j (A,), for all Ai £ L. (6) 

jeM 

For the sake of clarity, we may assume that L = M — {1, • • • TV) . He uses the notation : 


fi • • ' fk 

A\ • • ■ Ak 


f 


for denoting the coefficients of fk in the development (6). Then, in [12], he shows that if one 
has a Cebysev system (theorem 4.1 pp. 106) : 


f\ ' ' 

■fn 

f 


' fi 
A 2 * 

• fn- 1 

••A n 

/l- 

' fi’ 
Ax- 

• fn- 1 

■ A n _x 

/] 

A x - 

• • A n 


’ fi" 
_ * 

‘ fn- 1 

••A n 

/]- 

' fi- 
Ax- 

• • fn- 1 

• • A n _x 

/] 
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This expression is a direct generalization of the classical Newton formula. Let us make 
several comments on this formula when applied to our problem : 

1. If one adopts the lexicographic ordering 1,X,Y,X 2 ,XY,Y 2 ,- • • ,X n ,X Y, 
• • -,XF n_1 , Y n , the previous formula (7) must be applied n + 1 times to go from a 
total degree n to a total degre n + 1. One must also store quite a lot of terms like 


to build the divided difference table. For example, to go from degree one to degree 
two, one must evaluate and store approximations of gradients by means of approx- 
imations on triangles, combine them to obtain approximations on sets of four points 
(Cg sets), of five points (C| sets) and of six points (one set). Moreover, to go from 
approximation on k points to k + 1 points, (k + 1) x (k + 1) determinants must be 

evaluated. 

2. FYom a numerical point of view, the basis X'Y> or ( X - a)‘(Y - b)‘ are not well suited 
to calculations. This can easily be seen since for any pair («,&)> 



i + j <n 
(x t ,yi) € <5 (n) 


= [(zi - a) < ( yi - &) J ] 


i + j<n 

{xi, yi) € <S (n) 


If (a, b) is any point of S (n) and if h = max (xiiyi)eS w{\ x i~ a l, then Hadamard s 

inequality shows that : 

|A 5 (n)| <h^\ 

where k(u) = 1 + P (P + ^ + ^ = °(n 4 ) that one reaches very quickly 

machine zero though the linear system may be well conditioned. 


An alternative to this last point is to use local coordinates such as the barycentric ones. In the 
E.N.O. method we will develop in a further section, for each point, the natural barycentric 
coordinates are not known a priori so that the work has to be repeated at each interpolation 
call. If this is included in an iterative algorithm, the cost (and the storage) seems to be much 
too important at least for the cases we have considered in this report. For all these reasons, 
we have preferred to use classical inversion techniques for linear systems. 


2.2 A recurrence formula 

In this section, we wish to show another recurrence formula that enables us to obtain the 
coefficients of the expansion of total degree n from those of the expansion of total degre less 
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than n in only one step. This recurrence formula may be viewed as another version of that 
given in [12], theorem 3.1, page 400 and will be useful in section 2.3. 

Let us begin with some notations. Let P u •••, P N{n) be a basis of IR^F] such that 
Rit ' ' 'i p N( P ), p < n, is a basis of IR^ [A" , T] . The three bases we have considered in section 
2.1.2 are of that kind. Let Ai, • • •, be admissible points. We set : 

Rl=(P: i(A L )---P N{n) (A L )) T 

so that the solution of the Lagrange problem (where the u( s are given), 

«. = Yh a j p j(Ai) for all Ai € L, 

l<j<N(n) 

may be seen as the solution of the linear system 

M (oti- • ■ ajv(„)) = U = (ui ■ • • ujv( n )) , (8) 

where the Lth row of A4 is Ri,. The solutions of system (8) are : 

n detjRl • • • - ^AT(n)) 

L det(R 1 ...R L ...R N{n) y ^ 

where Rl = U denotes the L-th column. 


Lemma 2.1 Let (A)i<»<W(n) be an admissible set in which any of its N(p),p < n subsets 
is admissible. Then, for a given p < n, let I = {z’i, • ■ • *';v(p)} and J be ordered sets such that 
I\JJ — {!,•■• ^(n)}. If A = (a,- j ) is a N(n ) x N(n ) matrix, we set 


and 


det(A)i = det(a , j) 


1 < i < N(p) ’ 


det(A)j = det(ai j ) 


| j el 

N(p) + l<i< N(n ) ' 

j e J 


Let us also denote by ctyj the coefficient of Py of the Lagrange problem of degree p for nodes 
in I (for degree n, we omit the subscript I ). 

Then for any L' < N(p), we have 


A n ) _ v \L L' Av) 

a L ~ l^I,card(I)=N(p) A / a L' I 

\L V _ (_f \<r(l) det(Ri • • • Ry ■ • • Rfr r/)/ det(i?7v(p)+i • • • Ry • • • Rn)j ^ ^ 

1 K } det{Ri • ■ • R n ) ’ 

where a(I ) = 1 + p(p + l)/2 + & n (10), Ry appears a first time at the L'-th row of 

det(Ri • • • Ry • ■ ■ R n ,)i and a second time at the L-N(p)-th row o/det(f?AT(p )+ i • • • Ry • • • Rn)j- 
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Proof : By switching columns L and U in (9), one gets : 

( n ) _ det(Ri • • • • • • Ry • • • ^yv(n)) 
aL det(Ri • • • Rl • • • Rn{u)) 

Then, a direct application of the generalized Lagrange formula (see [13], pp. 19-22) to the 
previous expression gives 

det(R \ , • • • Rl * * • Rl< • • ■ Rn(h)) = 

E/,cord(/)=7V(p)(- 1 ) <T(/) • • • Rl“ • Rn(p))i det(R N (p)+i Rv Rn{v))j 

Since any iV(p) subset is admissible, one has 

( p ) _ detfRi ■ • • ‘ ' * ^7V( P ))j 

ai ' 1 det(Ri ■ • • Rv • • ■ Rn(ti))i ’ 

and the results follows immediately. • 

Lemma 2.2 VFif/i the assumptions of lemma 2.1, then, for p < n, if N(p ) < L < N(n ) and 

V < N( P ), 

E A / L ' = o. 

/,card(/)=N(p) 

Proof : Apply the Lagrange formula to 

det{R\ • • • i?L' • ■ • Rn(p)Rn(p)+i • • * Rl 1 Rn(h ) ) = 0> 

and interpret the coefficients in terms of a’s, all equal to 1, and A’s. The lemma 2.1 gives 
the result • 

2.3 Approximation of smooth and unsmooth function by poly- 
nomials 

The problem of interest in this section is the following : Let u be a real function defined 
on an open subset of 1R 2 . We assume that u is n times continuously differentiable on fl 
except perhaps on a subset of ft consisting of a finite collection of locally C 1 curves. Let now 
T be a mesh. For each point of T, we consider a Lagrange interpolation of u. Is it possible 
to localize the regions of smoothness of u from the coefficients of the Lagrange interpolation 
of u ? The answer is yes, at least for second and third order approximations if additional 
assumptions are made on the mesh. These assumptions guarantee that one can solve the 
Lagrange problem for any order from 1 to n, and may be seen as a very natural generalization 
of classical conditions used in the finite elements theory [14], 

For functions defined on 1R, one knows that the divided differences of u satisfy : 
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• If u is smooth on an interval I containing x k , ■ • • x n , then there exists (El such that 

[x I ,x s ,---,x n | U ] = 

n\ 

• if has a jump [u^] on 7, one has 

[xi,x 2 , ••■,*„!«] = 0([u w ](x„ - Xi)^ — n ). 

In this section, we intend to generalize these relations, and in particular, the second one 
since this problem seems (surprisingly) not to have been studied yet. The proof appears to 
be technical, and we have not been able to prove it for any total degree. The proof is divided 
into two parts. In the first part, we study the case of a stencil S M of N(n) points where u 
admits two values, 0 and 1. We show here that the polynomial of degree n that interpolates 
u is exactly of total degree n. Then, we define a condition on the stencils that appears to be 
a generalization of the one that says that triangles must not have too small angles to ensure 
a uniform error bound for classical finite elements [14]. Then, using Lemmas 2.1 and 2.2, 
we obtain our result. Let us begin with the case of a stencil in which the convex hull u is 
smooth. 

2.3.1 Case of a “smooth” stencil 

This problem has been studied by for example Ciarlet and Raviart in [15]. Let us recall one 
of their main results : 

Theorem 2.3 Let E = {a,}^ be an admissible (for degree k) set of points oflR n , and let 
h and p be respectively the diameter of £ and the suprenum of the diameters of the spheres 
contained in the convex envelop K o/£* Let u be a function that admits everywhere in K a 
k + Ith derivative D k+l u with 

M k+ i = sup{||Z) fc+1 u(x)||; x E I<} < +oo. 

If P is the unique interpolating polynomial of degree < k of u, we have for any integer m 
with 0 < m < k, 

h k+1 

sup{||D m u(x) - D m P(x)\\-,x E K } < CM k+l — , 

P 

for some constants 

C = C(n, k, m, £). 

Moreover, if S' is obtained from S by an affine transformation, then C(n,k,m, E) = 

C(n, k,m, £'). 

From this inequality, one sees that the “flatter” E is, the poorer the estimation is. A 
direct application of this theorem gives a generalization of our first statement. 
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2.3.2 Case of an “unsmooth” stencil 


Study of a simplified problem 

Let us consider S^ an admissible stencil of cardinality j -l a nd So, Si two non- 

empty subsets of <SbO having empty intersection, the union of which is S^ n f Let us consider 
a polynomial P of total degree n such that for all points of So, P has value 0 and for those 
of S\, P has value 1. Then we conjecture that : 

Conjecture 2.1 If is admissible, then the total degree of P is exactly n. 

We have not found in the literature any general proof of the statement, but we have the 
following lemma : 

Lemma 2.4 If any subset of cardinality ^ , k < n, of is admissible, then 

the conjecture 2.1 is true for n = 1,2 

Proof : The proof will be given for degree 1,2. We also indicate the difficulties for higher 
degree. 

• Degree one. The stencil is made of three points A,B,C that form a triangle. P is 
either of type or 1 — and is of degree exactly one. 

• Degree > 1. Let us assume that P is at most of degree n — 1. Set N = card(So) and 
M = S\. We have N + M = N(n). One may assume that N < M by changing P into 
1 — P. So, 2N < N(n). In the following table, we give the maximum values of N up 
to degree 6 : 


degree 

(ft + 1) (ft + 2) 
2 

N 

iy max 

2 

6 

3 

3 

10 

5 

4 

15 

7-8 

5 

21 

10-11 

6 

28 

14 


From this table, one can see that there is always, for degree 2, at least 3 points that 
have the same value. Since these three points are admissible for degree one and since 
by assumption, P is either a constant or of degree 1, we see that it must be a constant 
which is absurd since it takes two different values. The same argument applied to 
degrees n=3,4,5 shows that P must be of degree exactly n — 1 because we always have 

(n + 1) (n + 2) ^ n (n + 1) 

4 2 
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but fails for degrees greater or equal to 6 (because the previous inequality does not 
hold if n > 5).« 

With this in hand, we get the following result, if Ri j is the following vector : 

Ri j = ((xi - x 0 )'(j/i - yo) ] • • • (xn - xo)*(yN - yo) 3 ) ■ 

Lemma 2.5 Let (x 0 ,yo) be any point of the convex hull E of S^ n K Assume that conjecture 
2 A is true . Let P be a polynomial that is 1 on S\ and 0 on So where both sets satisfy the 
assumptions of lemma 2-4- If h — max{(\xi — xo\,\yi — yo|) 5 (^/, Vi) £ S^}, and if S ^ 
satisfies, for some a > 0, 

Mm{xo6S,Me([ 11 ^j 1 -i 1 ^ji]|}>a, (11) 

then there exist two constants C\{n,ot) and C^fa) such that the coefficients of the Taylor 
expansion of P around (zo,J/o) : 

t+j<n 

satisfy 

C 2 (n)h~ n > £ l a » j\ Ci(n,a)h~ n . (12) 

i+j=n 

Remark : 

L The set of points we define is clearly not empty because, on the convex hull of S ^ 
(which is compact), the function defined by the right hand side of (11) is continuous, 

2. If <S (n) were a triangle, for example, the minimum of that function would be the mini- 
mum of the absolute value of the sines of its angles. 

Proof : We adopt the lexicographic ordering of monomials. The set of admissible stencils 
satisfying condition (11) and |ft| < C is a compact subset C of ]R;' v '( n ). Let us consider the 
real functions defined on C by : 


<f>i j{A\, • • • > An) — X) \det(Roo -Ui f?on)| > 

t=l T n 
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Let also V, j be J2i+j<n P(Aij)Uij where the “1” is at the L-th position (we refer to the 
lexicographic order). Let us note that P(Aij) is zero or one, and its value depends only on 
ij and not S^ n \ 

It is clear that iftij is the coefficient of ( X — Xo)’(K — yoY in the Taylor expansion of P. 
It is also clear that 

h l+J \det[Ro o • • • Ri j • • • Ro n ]| < <t>i j, 
by the triangle inequality. So, 


h n 6(Ai,- • - ,A n ) = h n £ ki| > E 

*+i=n t+j=n 


detjRoQ • • • , Kj • • • Ron) 
fit j 


The left hand side of this inequality is a continuous function on C and hence reaches its 
minimum. This minimum cannot be zero because that would mean that all of the a; fs are 
zero which is absurd by assumption, and from lemma 2.2. 

Now let us turn to the second inequality. We have 


h'e{A u ---,A N ) = h~ £ !«,,!= £ |fcj|. 

i+j=n i+j=n 

The left hand side is also a continuous function on C and is bounded above. Clearly, the 
latter constant does not depend on a. • 

Corollar 2.6 With the assumptions of lemmas 2.4 and 2 . 5 , let n and p be integers satisfying 
p < n. Choose L and L 1 such that N(p — 1) < V < N(p ) and N(n — 1) < L < N(n). Then 
there exist two constants C\ and C 2 such that 

C 2 {n)h~ n+P >X^ L ’> C 1 {n,a)h~ n+P , 


for any subset of cardinality p. 

Proof : We apply the definition of Af L ' and the same techniques used in the previous 
lemma • 

Now, we may state our main result : 

Theorem 2.7 Let be a stencil satisfying the assumption of Lemmas 2-4 and 2.5 for 
n = 1,2. Let u a real function defined on an open subset ofLl in 1R 2 being C 1 except perhaps 
on a locally C 1 curve where its n-th order derivative may have a jump [D n u] such that the 
intersection T of that curve and the convex hull of S ^ is not empty. Then there exits a 
constant C(n,a ) > 0 such that the coefficients in the Taylor expansion towards a point of T 
satisfy : 

E k j\ > C(n,a) 1^1. (13) 

i+j=n n 
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Proof : Let us assume that u admits p-th continuous derivatives but its p + 1-th ones have 
a jump on a locally C 1 curve C. Let be the coefficients of a Lagrange interpolation of 
degree p on a suitable subset of The recurrence formula of Lemma 2.1 gives : 

«5 l = E 

I,cardI=N' 

Let S+ and <5_ be the subset of S ^ defined by 


a £ 5 + ( resp. S- ) iff D{ju(a ) — ► Dfj( resp. D tJ ) 


These two sets exist because of the topological nature of the curve C (see Figure 1). 

Let t > 0 and x 0 be a point of J. There exists tj > 0 such that for max{\a x — xo|> \a x — 
x 0 |} < 77 , we have the following development of u for points of <S + : 

u(a) = E E Diju(x 0 )(a x - x 0 ) i (a s/ - y 0 ) j + E ( D t> + ‘K 1 )) ( a * _ x ° Y( a v ~ 

/=0 i+j=l ? +i == p + 1 


where a x and a y are the x- and y- components of a and Df is the limit of the i j partial 
derivative of u when x — > x 0 on the right and |o(l)| < t . The same type of development is 
also true for points of S _ ; Df is replaced by D~. The properties of the determinants allow 
us to get for any subset of cardinality N(p ) that, by changing the ordering of the points if 
necessary, 


Df{al - x 0 )*’(«l - Vo) 3 
Dtj« ~ x oY( a y -yo) j 

Rqo ••• ~*o )’« +1 -yo ) 3 ■■■ Ro V 

D-Aaz {p) - x oYK {p) -yo) 3 

det{Rm • • • Rof) 

o(l)(<4 - x 0 )‘(aj - yo) 3 

o(i)( a ™ - zo) l « - yoY 

Roo ••• °(l)( a ™ +1 -x o y{a™ +1 -y 0 ) j i?o P 

o(f)(a^ (p) - xp )*(af (p) ~ VoY 

det(Roo ■ ■ ■ Rq p ) 


(14) 


where the index V stands for (z, j) in the lexicographic order. For the sake of simplicity, let 
us denote by pi and v\ the coefficients obtained from the first part of the right hand side of 
equation 14 by replacing Df by 1 and D~ by 0 for pi and vice versa for vi so that 

«L'J = [Df + 0(1)) fl! + (Df + 0(1)) Vj- 


12 


We have to notice that the sum of /z/ and Vi is one, and that if all the points of / are on the 
same side of C, then either /z/ or zq is zero. Then, using lemma 2.2 and the above remarks, 
we have, if [Ay] = A+ — A~- : 

QL = [Ay] X L + X L + X i 

card(I)=p card(I)=p card(I)=p 

so that 

|[z><j]| E + £ E (M + M)l*“'l- 

card(I)-p card{l)-p 

We can also get a similar inequality by replacing /z/ by vj if the points of I are not on the 
same side of C so that : 

h P ~ n I [Ay] | Ecard(/)=p L < h p n {|[Aj]| \j2ca.Td(I)=p | 

+ |[Aj]| |Ecarrf(7)=p V I^I L |} 

+ 2 h p - n \a L \ + 2 eh p ~ n £ cord( /)= p (Im/| + M) l A / L '\- 

(15) 

When all the points of / are on the same side of C the factor 2 is replaced by 1. Because of 
inequality (11) and from Hadamard’s inequality, one has : 

,/EK- *o)>i - ^ + ,1 E («i - *o)“K - w) 2i 

V a 5+ \ a S- 

M + N * — ® 

This latter expression is bounded above by a constant C because L corresponds to {i,j) 
in the lexicographic order, without any additional conditions on the geometry. So, the 
inequalities of (15), with the help of the first inequality of lemma 2.4 become (up to a factor 
2 sometimes) : 

| [Ay] | X A / V < ^ M + ^ 

card(I)—p 

for a given constant C' . Now, summing up these inequalities for i + j = n, with the help of 
the second inequality of lemma 2.4 one gets our result for e small enough. • 

This result enables us to detect the regions of smoothness from those where a jump in 
one of the derivatives occurs. It is true when conjecture 2.1 is true, and at least for degree 
1 and 2. 
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3 An E.N.O. Reconstruction Technique 

In recent papers, Harten and several other authors [3, 4, 5, 6, 7] have tried to derive numerical 
methods that are able to achieve a higher order of accuracy than classical TVD methods. 
There are several versions of these techniques, but they can be generally viewed in the 
following way : starting from some approximation of a real function u (point values or 
average values in some control volumes), find a pointwise high order approximation u. Two 
tools are then used : 

• an essentially non oscillatory Lagrange type interpolation of a function w, 

• this function w may be u itself if one starts from point values or either the primitive 
function of u or its convolution product with the characteristic function of a copy of 
the control volume if one can pass from one to another by a translation. 

The latter deconvolution technique can only be applied, at least in its standard version, to 
regular meshes as shown in section 3.2. 

In this section, we want to adapt both tools to the situation of unstructured meshes. At 
least for the second point, the situation seems at first glance very bad : the reconstruction 
via primitive function cannot be applied in the case of unstructured meshes because the only 
solution would be to apply it to integrals over domains like T>m = [ai, &i] x [a 2 , b?], possibly 
with a suitable transformation of the plane as in [8], which are not in general the union of 
control volumes. 

Now, the reconstruction via deconvolution technique can only be applied to regular 
meshes (i.e. meshes where the control volumes are translated form one node to another). In 
this section, we show how to adapt this technique for irregular meshes. 

In what follows, T is a triangulation of 0, a domain of 1R 2 , u is a function defined on 
that domain. Around each node i of T, we may define a control volume C, in many different 
ways. An example is (see Figure 2) the control volume whose boundary is the segment 
joining the centroids G of the triangles ( i,j,k ) having i as a vertex and the middles 7 a , J 2 
of the segments of those triangles (type I). Another type of control volume is obtained by 
considering each triangle as the control volume of its centroid. The triangulation we need 
to describe the ENO algorithm is not T but another one built from the centroids of the 
triangles of T . These control volumes always satisfy : 

• UigT C* is the numerical domain fi/,, an approximation of the physical one, 

• C, f| Cj is of empty interior when i / j (generally speaking, a collection of segments or 
an empty set). 
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3.1 Non oscillatory interpolation 

Let us consider n > 0. In this section, we show how to generalize the n + 1-th order E.N.O. 
interpolation technique exposed in [4, 5] to unstructured meshes. As the results of chapter 
2 indicate, we must deal with meshes where the stencils we need satisfy the property of 
lemmas 2.1 and 2.2. This will be the case for most meshes. 

Let T be a triangulation of ft, a domain in IR 2 , and u a function defined on that domain. 
The results we have obtained in chapter 2 can be summarized as follows : 

• if S W is an admissible stencil such that u is smooth in its convex hull, then 

Y l a * il> 

i+j=zn 

remains finite, 

t if in the convex hull of S^ n \ u admits continuously differentiable derivatives only up 
to the order k < n, then 

£ k si = o« u <">] 

t+j=n 

where h is the diameter of S^ n \ 

Then, as suggested by [4, 5], the E.N.O. algorithm we propose is to consider IIi(u) defined 
on Ci by the following recursive algorithm : 

For a node i , 

1. Let {T;} be the set of triangles of T having t as a vertex. Consider all the linear 
interpolations where the T,’s are the stencils. Choose the one, f mm , where the sum 

Y K ;l’ 

i+j - 1 

is minimal. We set <S^ = the nodes of T m in- 

2. Let 5O-1) be the stencils defined at the previous step. Consider all the nodes sur- 

rounding 5O- 1 ) in T and consider all the stencils obtained form S (n-1) by adding n + 1 
of the nodes surrounding Choose the stencil minimizing : 

Y 

t+j=n 

We have intentionally left the second point imprecise because it is obvious that the 
number of stencils to consider is in general huge. To give an example, if one considers Figure 
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3 for which n — 2, one sees that possible stencils are the vertices of triangles T m i n and 3 of 
the ten other triangles. This can be repeated for each of the three edges of T mtn and leads 
to a total of 3 x Cf 0 = 360 possible stencils ! So, one has to define criteria for choosing the 
“good” and “bad” stencils. These criteria are essentially heuristic and a priori ones . 

One that seems natural is that when one considers the control volume around each node, 
the collection of the control volumes of all points of the stencil should be convex. Another 
one is that the criteria leads to the smallest possible number of stencils, but the stencils must 
not be confined in a particular angular area of the plane, in order not to favor any direction. 
With this in mind, two possible sets of stencils for third order interpolation are: 

• the nodes of triangle T mi > „ plus, for each of its edges, the three additional nodes of 
triangles 7\, T 2 , T 3 . This leads to a maximum of three stencils per triangle, 

• or the nodes of triangles of T mtn plus, for each of its edges, the three additional nodes 
of triangles of 

- T x ,T 2 ,T 3 , 

- Ti,T 2 and T 4 or T 5 , 

- T x ,T 3 and T 6 or T 7 , 

- Ti, Tg or Tio and one of the six triangles T 2 , T 3 , T 4 , T 5 , T 6 , TV. 

The second solution leads to a maximum number of 52 stencils once T mtn has been found. 
We have made several tests to evaluate the “performance” of each type of stencil. They are 
given in section 4. 

This particular interpolation is n + 1-th order accurate; because u is a polynomial of 
degree less than n, we have fli(u) = u. This property ensures the n + 1-th order accuracy 
[15], and in particular, we have the estimations of theorem 2.3. 

3.2 Deconvolution technique revisited 

If is a regular mesh of R and u is a real valued function on IR, the reconstruction by 

the deconvolution technique consists of applying the previous algorithm, not to u but to 



where, as usual, the mesh size Ax = Xi + i — Xi is constant and x; +1 / 2 = x; + Ax/2. In 
particular, we see that u,- does not depend on i and that u(x;) is the average of u on 
l x i-i/2>%i+i/2]- These values are assumed to be known. Let IIi(u) be the m + 1-th order 
Lagrange interpolation as decribed in the previous section, with m > n. Then, the idea is to 
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perform a Taylor expansion of u and its successive derivatives around x,, to truncate them 
at order n — k, to replace the values u , u ^ by IIi(u) and its n successive derivatives, 

and to replace the values of u by those of Il 2 (u), the approximation of u we are looking for : 


n i(u)(s,-) = ep=i q|— 


= n=i k 


ni(u) (n) (x,) 



(17) 


where 

r i,2 (x- Xi y dx 

l ■'*<- 1/2 

1 " Ax 

The linear system is easily invertible because the matrix is upper triangular and its diagonal 
consists of all “l”s. Furthermore, it is shown in [4], for example, that the average value of 
Il 2 (u) over [x,_i/ 2 , x, + x/ 2 ] is exactly u,. Last, this approximation has the desired order of 
accuracy when u is smooth because polynomials are left invariant by the construction. 

This latter point is the fundamental reason why one achieves the expected order of 
accuracy . Polynomials are left invariant by the construction because the shape of the 
control cells does not change from one point to another. If this where not the case, i.e if 
Axj_j_i /2 = x t+ i — x, were not constant, the formula (16) would indeed depend on the point 
x,- and this property would be lost. In order to show this, we simply consider u(x) = x and 
a m + 1-th order interpolation that has values u,- at points x,-. Assuming that {xo, Xi, • • •} is 
the stencil selected by the E.N.O. algorithm, we have : 

IIi(u) = u 0 + K(x - x 0 ) -| where K = j- + ^zL 1 1 _ 

1 Z\Xj/2 


When the mesh is not regular, K ^ 1 in general. To obtain Il 2 (u) = u, one must have 


n i(w)(i/) = u(y) + aiu'(y), 
n 2 («)(y) = I< = u'(y). 

The second equation indicates that one must have K = 1 which is, in general, not true. 

To overcome this problem, we propose the following technique : apply the ENO search 
algorithm not to u defined by equation 16 but to : 

u(y) = 777 u(x + y-x 0 )dx, (18) 

area(C < 5(n)J 7c s(n) 
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where 5 W is any possible stencil around the node x 0 that one has to test and C S { n) is the 
union of the control volumes of each node in S(«) : 

<7*0= U Cj. (19) 

j€S( n ) 

Now, one has to evaluate the integral (18) from the average value of u. This cannot generally 
be achieved for any function, but is possible for the polynomials of IFUipf, Y], at least in 
general. This will be true for the N(n ) linear forms over IRnfX, Y\ because 

<P>i= — [ P(x)dx, forall PGlRn[X,F], (20) 

areaGi JCi 

are independent. For all the meshes we have considered, these linear forms where always 
independent, so the problem had a solution. If this is true, then one can find coefficients 
a i{y)i 1 < / < N(n) so that 

1 r N{n) 

— r / u(x + y- x 0 ) dy='£ / a t (y ) < u >,, (21) 

area{C^(n)) dc s ^ 

when u belongs to IR ^[X,Y]. If not, the equation (21) gives a n + 1-th order quadrature 
formula. With all this, we get the following theorem whose proof is obvious. 

Theorem 3*1 The algorithm defined by the E.N.O. technique with (18), (19), (20) and (21) 
leaves invariant the polynomials o/IRnfX, Y] and hence gives an + 1 -th order approximation 
of smooth functions . Moreover , this approximation is conservative up to the quadrature 
formula 2L 

3.3 Some remarks for the practical calculation of the reconstruc- 
tion 

In section 2.1.3, we have discussed the problem of the practical determination of the coeffi- 
cients of a Lagrange interpolant because the linear systems to be solved have in general small 
coefficients. The same problem also arises here if one uses the monomials (X — xo)*(K— yo) 3 to 
determine the ai(y ) of equation 21. This can easily be seen by using Hadamard’s inequality 
as in section 2.1.3. 

For a given node xo, the E.N.O. technique we propose naturally introduces one triangle 
having that node as vertex, the triangle T m ,„ as in Figure 3 . So, as in section 2.1.3, we will 
use the barycentric coordinates towards that triangle for practical computations. 
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4 Numerical examples 


We have performed several tests on the second and third order E.N.O. interpolation and 
E.N.O. reconstruction, but we only report the third order results since they are a priori 
more challenging. In particular, we intend to check numerically that the expected order of 
accuracy is in fact reached for smooth functions. 

The two types of stencils have been tested on the latter case. The use of the second type 
of possible stencils results in a much more expensive approximation (in general, one must 
test 52 stencils per triangle versus only 3 in the simplest version) and the results have never 
been dramatically improved. All the results that are presented bellow have been obtained 
with the 3 stencil version of the method. 

In all these examples, we have assumed that the control volumes are of “type I”. The 
practical calculations of the averages in these control volumes have been performed with a 
5-th order quadrature formula [14]. 

The tests on smooth functions will be performed on : 

u(x , y) = cos(2w(x 2 + y 2 )). 

We have displayed in Figure 4 the L°° error of the interpolation. The plain curve with 
squares is obtained with the E.N.O. interpolation, the plain curve with circles is obtained 
with the E.N.O. reconstruction. The dashed line indicates the slope —3. One can see that 
the expected order of accuracy is indeed reached. The mesh size h has been measured by 
choosing the largest segment of the triangulation. All the error estimates have been obtained 
on irregular meshes as the one presented on Figure 5. These meshes are obtained as random 
perturbations of regular structured meshes. The set of points that one obtains is triangulated 
by the Brower algorithm to get a Delaunay triangulation. The main difference between such 
a mesh and the regular structured one is that the number of triangles each node belongs to 
is different. We also have done the same tests with regular meshes, and we have not seen 
any degradation of the convergence. 

The locally smooth function we have chosen is obtained by a modification of that used 
by Harten in [7] for example : if <j> is any angle, let be : 

f if r < - 3 , U{x,y) = -rsin(|r 2 ), 

/*(*, y) = < if r > §, U(x, y) = 2 r - | sin ( 3 ?rr), where r = a; + tan (<f>)y, ( 22 ) 

l if M < §, U{x,y) = |sin(27rr)|, 

and let u be : 

if x<^ cos ny, u{x,y) = fr^(x,y), 

Si * 

if x >2 cos Try, u(x, y) = f_^p{x, y) + cos {2-iry). 
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The function defined by (22)-(23) shows discontinuities in the function itself and its first order 
derivatives; some of the discontinuities are straight lines (never aligned to the mesh), one is 
a curved line where the jump changes from one point to another. Last, the behaviour of u is 
basically one-dimensional on the left of the curve x = costtj//2 and really two-dimensional 
on the right. 

A plot of this function is given in Figure 6. One should obtain straight lines and smooth 
discontinuity transitions contrary to what is shown in the Figure : this is an effect of the 
graphic device adapted to P\ interpolation. 

In Figures 7 and 8, we have displayed the node values of the E.N.O. reconstruction for 
two meshes (1600 nodes and 6400 nodes). To better see the behaviour of both approximation 
techniques, we also present cross-section on three lines : Y = 0.75, Y = 0 and Y = —0.45. 
The approximations are obtained from the 1600 nodes mesh (see Figure 5). The latter 
line goes through one of the triple points (see Figure 6). One can see that the various 
discontinuities and the smooth regions are well captured by both techniques. 

To end this section, we must note that the algorithm for choosing the stencil may lead 
to some difficulties at the boundaries as can be seen in Figure 6 on the left upper corner : 
the most left upper triangle of the mesh (Figure 5) does not admit any additional points of 
the type we consider to make a stencil. 

5 Conclusion 

In this report, we have developed two methods for the reconstruction of a function admitting 
discontinuities only on regular planar curves from their node values or from their average in 
control volumes that surround them. In order to give a firm basis to the Essentially Non- 
Oscillatory interpolation, we have studied the behaviour of the highest order coefficients of 
the Lagrange interpolation of smooth functions and unsmooth ones for which the discontinu- 
ities lie on regular curves. We have also given an adaptation of the so called “reconstruction 
via deconvolution” to irregular triangulated meshes. 

These techniques have been shown to work quite well on smooth and unsmooth functions. 
In particular, we have shown in these examples that the minimum number of possible stencils 
was sufficient for our purpose. 
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Figure 4: L°° error for f(x,y) = cos [27r(x 2 + y 2 )]. Squares : E.N.O. interpolation only, 
Circles : E.N.O. + reconstruction. Dashed line : slope -3 
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Figure 5: Typical mesh. 1600 nodes, 3042 triangles. 
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Figure 6: Exact function, Mesh with 6400 nodes. Min=-1.331, Max=2.650, 8 = 8.292 10 -2 . 
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Figure 7: E.N.O. reconstruction with the 1600 nodes mesh. Min=-1.308, Max-2.651, <5 - 
8.249 10" 2 . 
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Figure 8: E.N.O. reconstruction with the 6400 nodes mesh. Min— -1.325, Max— 2.650, 8 — 
0.8281 lO' 2 . 
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